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PREFACE 


In  sonar  evaluation  studies  there  are  at  present  three  approximate 
methods  of  determining  sound  propagation:  ray  theory,  mode  theory  and 
semi-empirical  formulation.  Under  the  usual  environmental  conditions, 
one  or  a  combination  of  two  of  these  methods  will  produce  satisfactory 
prediction  of  intensity  versus  range.  But  under  extreme  environmental 
circumstances,  favorable  or  unfavorable,  range  prediction  is  not  re¬ 
liable  with  any  combination  of  these  methods.  Another  limitation  of 
the  existing  methods  arises  with  the  increasing  reliance  on  signal  pro¬ 
cessing.  Static  methods  of  range  prediction  are  becoming  less  useful 
even  under  favorable  circumstances. 

This  Memorandum,  part  of  a  larger  study  conducted  for  the  Advanced 
Research  Projects  Agency,  presents  a  technique  which  may  prove  useful 
under  difficult  environmental  conditions  and  in  situations  in  which  a 
dynamic  solution  is  required.  This  latter  capability  should  be  of 
special  interest  to  those  who  are  concerned  with  sonar  development  in 
particular,  and  with  acoustic  propagation  in  general. 
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SUMMARY 


A  computer  solution  of  the  wave-dif ference  equations  Is  found  by 
using  an  interlacing  explicit-implicit  scheme.  In  the  computation,  the 
entire  two-dimensional  field  is  found  as  a  function  of  time. 

The  examples  considered  irvolve  propagation  in  a  homogeneous 
shallow-water  channel  where  tho  effect  of  superposition  of  discrete 
spectra  produces  characteristic  modal  patterns.  The  spatial  sound 
pressure  fluctuations  are  represented  as  plot -density  variations  along 
the  two-dimensional  channel.  Examples  of  discrete  propagated  modes  as 
well  as  evanescent  modes  are  presented.  The  size  of  the  field  that  can 
be  presented  is  limited  by  the  size,  accuracy,  and  speed  of  the  com¬ 
puter. 
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I .  INTRODUCTION 

In  this  study  the  difference  equations  representing  wave 
motion  are  solved  by  an  interlacing  explicit-implicit  procedure. 
This  scheme  in  its  present  form  was  applied  with  remarkable  suc¬ 
cess  by  J.  J.  Leendertse  to  the  solution  of  long-r  riod  ocean  wave 
motion/  *  Reference  1  demonstrated  that  this  method  is  capable 
of  producing  accurate  stable  solutions  of  the  nonlinear  wave 
equation  without  the  usual  restricting  assumptions.  The  present 

study  is  concerned  with  the  two-dimensional  linearized  undamped 
wave  equation. 
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II.  DIFFERENCE  * '"ATIONS 

The  basic  first-order  undamped  differential  equations  which  de- 

(2) 

scribe  the  seur.d  field  are 


|W  +  j»(c_SJ.  .  o 

p  St  Sz 


as.  +  e  as  +  p  as  .  o 

St  Sx  p  '.2 


where  the  quantities  §,  U,  W  represent  the  linearized  departures  in  the 
variables,  density,  and  two  components  of  particle  velocity,  respectively. 
In  this  Memorandum,  the  sound  velocity  C  is  considered  only  as  a  function 
of  the  depth  z;  the  density  p  is  a  constant.  Both  of  these  quantities 
could  be  functions  of  spatial  as  well  as  of  time  coordinates. 

The  difference  equations  which  represent  the  differential  Eqs .  (1), 
(2)  and  (3)  could  be  formed  in  a  variety  of  ways.  The  scheme  chosen 
here  uses  a  two -inter lacing  grid  system  in  which  the  density  point  is 
surrounded  by  two  stream  points,  as  indicated  in  Fig.  1.  Position  in 
th  grid,  formed  by  thsse  points,  is  indicated  by  two  subscripts.  The 
first  refers  to  the  row  and  the  second  to  the  column. 


U '  -  u  -  -I-  c  ^  ('■  ‘  -  ~  1  ) 

n,m  n,m  l p  n  °n,rrrfl  Jn,nr 


/  T  2  2 

W  =  W  ■  7  (C  ,.  5  ,.  -  C  l  ) 

n,m  n,m  £p  '  n+1  n+l,m  n  n,m 


(5) 


Physical  boundary,  free  surface 


Fig.  1— Scheme  for  the  boundaries  and  the  grid  location 
of  excess  velocity  and  density  in  the  channel 
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il  =  §  -  (u;  -  u'  .)  -  —  (w  -  w  ) 

:,m  n,m  l  n,m  n,m-l  l  n,m  n-l,m 


(6) 


The  primed  superscripts  designate  the  quantity  at  t  +  t. 


arid 


Note  another  equally  valid  formulation: 

U'  =  U  -  T~  C2  (S  . ,  -  ?  ) 

n,m  n,m  £p  n  n,m+l  n,m 

w'  =  w  -  (c2  § '  -  c2  ?'  ) 

n.m  n,m  j?,p  n+1  n+i,m  n  n,m 


§'  •  i  -  (U  -  U  .)  -  Si  (w#  -  w'  ,  \ 

n,m  n,m  l  n,m  n,m-l  x  n,m  n-l,m' 


This  formulation  is  used  in  the  second  half-time  step  wtien  the  alter¬ 
nate  solution  direction  is  used. 
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III.  THE  ALTERNATING  DIRECTION  IMPLICIT-EXPLICIT 
METHOD  OF  SOLUTION 

During  each  of  two  successive  time  increments  T,  a  different 
operation  is  performed.  Equations  (1)  and  (3)  are  first  solved  im¬ 
plicitly;  Eq.  (2)  is  then  solved  explicitly  in  the  first  time  inter¬ 
val.  In  the  second  time  interval,  Eqs.  (2)  and  (3)  are  solved  implic¬ 
itly  and  Eq.  (1)  explicitly.  These  two  operations  comprise  a  one-time 
step  of  2".  From  a  knowledge  of  the  boundary  conditions,  tt:o  two  groups 
of  implicit  equations  are  solved  by  eliminating  the  urk-  owns . 

The  alternating  direction  method  car.  be  applied  directly  to  Eqs. 

(4)  through  (6)  if  the  equations  are  put  in  the  following  form: 

First  half-time  step:  implicit  solution  over  the  row  n 


7PU'  .  +  §'  +jpu'  =  A 
"■  n,m-l  n,m  £  n,m  n,n 


B  =  U 
n,m  n,m 


First  half-time  step:  explicit  solution  over  the  column  m 


v'  =  W 

n,m  n,m 


(9) 
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Second  half-time  step:  implicit  solution  over  the  column  m 


•  7  p  H  .  +7  P  w  -a 

X  n-l,m  n,m  l  n,m  n,m 


(10) 
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Second  half-time  step:  explicit  solution  over  row  n 


u"  =  u' 

n,m  n,m 


(12) 


During  the  first  half-time  step  t  to  t  +  t,  Eqs .  (7)  and  (8)  are 
solved  implicitly,  while  Eq.  (9)  is  solved  explicitly.  The  implicit 
formulas  for  each  n  during  the  first  half-time  step  can  he  written  in 
general  terms 


-  P  U'  +  Q 
n,m  n,m  n  >  ni 


(13) 
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/ 

n  ,rr.-l 


-  R  ,  §'  +  S  . 

n ,m- 1  n ,m  n ,m-l 


(14) 


7 


where 
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7  p 


n,m 


+  i 


PR, 


h  ,m-l 


+  7  P  S_ 


n.m  -  .  t 

1  +  7  p  R  . 
I  n,m-l 
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JLC2 
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pi  n 


P 

n  ,m 
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S 

n  ,m 
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r.  ,m 


T 

Pf 


P 

n,m 


if  m  4  1 


if  m  -  1 


U/  represents  an  external  driving  function  at  the  zero  position;  it 
n ,  L 

is  zero  otherwise.  In  tlie  following  time  step  t  +  t  to  t  +  2t,  Eqs  .  (10) 
and  (11)  are  solved  implicitly,  while  Eq.  (12)  is  solved  explicitly. 

The  generalized  form  of  the  equations  is 


P  w  +  q_ 
n,m  n,m  Ti,m 


(15) 


t- ' 

n-l,m 


r  i  5  +  s  . 

n-i,m  n,m  n_i,m 


(16) 
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where 


p 

n,m 


**n,m 


r 

n,m 


s 

n,n 


Equations  (13)  through 
computation. 
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The  explicit  Eq.  (9)  is  then  solved  over  the  columns  m  beginning  with 
n  *  2. 
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In  the  following  half-time  step  p  ,  q  ,  etc.  are  evaluated  in 

n,m  n  ,m 

succession  over  each  column  beginning  with  n  *  2.  The  implicit  Eqs .  (15) 

and  (16)  for  §  ,  W  can  then  be  solved  in  the  same  way  as  the  first 

n ,m  n,o 

half-time  step,  except  that  the  boundary  conditions  now  depend  on  the 
surface  and  bottom  of  the  channel.  The  water-air  interface  cannot  sus¬ 
tain  any  change  of  pressure.  The  boundary  conditions  at  the  free  sur¬ 


face  of  the  channel  n  “  n  can  be  satisfied  if  £  ■  0.  Equation  (16) 

max  an  ,m  n  '  ' 

max 

becomes 


n  -l,r 
max 


n  -l,m 
max 


Beginning  with  this  equation,  at  the  top  of  the  channel  the  solution 

proceeds  to  the  bottom  of  the  channel  in  a  procedure  similar  to  that 

used  to  compute  the  first  half-time  step.  For  the  rigid  boundary  at  the 

bottom  of  the  channel  (Fig.  1),  the  velocity  in  the  vertical  direction 

is  W''  *  0.  Then  Eq.  (15)  can  be  written 

l  ,m 


^2,m  ^2,m  ^2,m  +  ^2,m 
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V.  DISCUSSION  OF  THE  EXAMPLES 


The  interlaced  implicit-explicit  method  was  programmed  in  FORTRAN 
for  RAND's  7040-7044  computer  for  a  rectangular  channel.  Various  grid 
sizes,  time  increments  and  frequencies  were  tried.  The  numerical  sta¬ 
bility  of  the  calculations  and  the  basis  for  grid  and  time  step  size 
have  been  discussed  in  Ref.  2.  For  the  examples  presented,  a  20  cps 
source  of  the  form 

A  sin  (iut  +  <p) 

was  applied  at  one  end  of  the  channel.  In  order  to  reduce  the  transient 
disturbances,  a  zero  phase  angle  was  used.  Peak  amplitudes  of  the  steady- 
state  solution  of  excess  density  at  each  data  point  of  the  channel  were 
retained  and  pi  >tted  on  the  S-C  4060  microfilm  recorder.  A  density  plot 
representative  of  the  peak  density  amplitude  was  constructed  by  making 
the  size  of  the  character  (X  in  this  case)  represent  a  density  amplitude 
range.  In  Figs.  2  and  3  the  darker  regions  represent  less  sonified  re- 

* 

gions,  while  the  white  regions  represent  those  which  are  highly  sonified. 
These  plots  represent  only  a  portion  of  the  steady-state  field. 

Examples  of  patterns  of  locked  and  evanescent  modes  are  shown  in 
Figs.  2  and  3.  As  indicated  in  Fig.  2(a),  the  source  was  placed  in  the 
center  of  the  end  of  the  channel.  This  resulted  in  the  energizing  of 
modes  1  and  3.  In  Fig.  2(b),  the  source  was  placed  near  a  nodal  point 
of  mode  3;  hence  mode  3  was  suppressed  to  some  extent  and  energy  was 


As  programmed,  the  maximum  peak  amplitude  is  divided  into 
five  equal  increments.  The  first  increment  (a  blank)  represents 
the  excess  density  in  the  range  1.0A  to  0.8A,  and  the  four  successive 
divisions  represent  decrements  in  amplitude  with  the  largest  X 
designating  the  range  0. 2A  to  0.0. 


Rigid  surface 


Fig.  2— Peak  amplitude  of  the  excess  density  in  a  channel 
237. 5  ft  deep  and  1550  ft  long 
(white  regions  are  highly  sonified) 
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directed  to  modes  1  and  2.  These  two  figures  represent  the  channel 
driven  at  a  frequency  below  mode  3  cutoff,  as  the  attenuation  in  the 
x-direction  demonstrates. 

In  Fig.  3  one  boundary  condition  -as  changed  to  represent  a  free 
surface,  as  described  in  the  preceding  section.  All  other  parameters 
are  the  same,  the  source  is  in  the  middle  of  the  channel,  and  modes  1 
and  2  are  propagated  while  mode  3  and  higher  modes  are  evanescent. 
These  results  are  very  much  like  those  of  the.  optical  fringe  analog 
of  sound  propagation  of  Ref.  3. 

It  is  possible  to  compare  these  results  with  ray  theory  by  veri¬ 
fying  that  the  horizontal  pattern  repetition  distance  A  given  by  ray 

(3) 

theory  is  the  same  as  the  A  indicated  in  Fig.  3.  From  ray  theory 


A 


(18) 


where  n  and  m  are  the  mode  numbers,  H  the  distance  between  two  density 
release  surfaces,  and  X  the  wavelength.  The  phase  change  at  the  rigid 
surface  is  equivalent  to  a  reflection  at  a  free  surface  lying  at  a  dis¬ 
tance 


H 


t 


=  _i_  X 

sin  0  4 


(19) 


below  the  rigid  surface.  Here  0  is  the  angle  that  the  ray  makes  with 
the  boundary  for  in-phase  regions  to  appear.  From  Fig.  3  this  angle 
is  approximately  20°.  The  enhanced  channel  depth  is  then 

H  “  H'  +  Hl 


(20) 
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From  the  parameters  assumed  in  Figs.  (2)  and  (3) 


sound  velocity  V 
wavelength  X 
channel  depth 
grid  size  4 
time  step  size  t 


=  5000  ft/sec 
=  250  ft 
=  237.5  ft 
=  12.5  ft 
-  0.00125  sec 


the  horizontal  repetition  length  calculated  from  Eqs .  (18),  (19)  and 
(20)  is  approximately  739  ft,  while  the  measured  distance  from  Fig.  3 
is  approximately  12.5  X  57.0  ft  *  712.5  ft.  The  difference  between 
the  ray  theory  calculation  and  the  calculation  in  this  work  probably 
can  be  attributed  to  quantization. 
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VI.  CONCLUSIONS 

It  is  important  to  recognize  the  distinctive  features  end  the  in¬ 
herent  limitations  of  the  implicit-explicit  method  compared  with  the 
methods  now  in  use.  It  is  difficult  to  infer  all  of  the  features  and 
limitations  of  this  method  from  these  simple  examples,  but  the  follow¬ 
ing  are  evident : 

FEATURES 

1.  The  entire  field  is  found  as  a  function  cf  time  in  the 
calculation . 

2.  Other  than  the  grid  size,  no  restrictions  are  placed 
on  the  velocity  field  and  boundary  conditions. 

3.  The  parameters  describing  the  field  can  vary  in  time. 

LIMITATIONS 

1.  The  field  size  is  limited  by  the  accuracy  and  the 
storage  capacity  of  the  computer. 

2.  The  computation  time  is  a  function  of  the  field  size 
and  the  computer  fast-storage  capacity.  For  example, 
one  hour  of  RAND  7 O'* 0-7044  computer  time  and  most  of  the 
fast-storage  capacity  were  consumed  in  pi.djcing  each 
solution  as  shown  in  Figs.  2  and  3. 

It  is  generally  accepted  that  ray  tracing,  normal  mode,  or  semi- 
empirical  oKthod,  or  some  combination  of  them,  is  adequate  for  most 
naval  applications,  orovided  that  the  environmental  parameters  and 
sonar  system  parameters  are  available  and  are  used  correctly.  How¬ 
ever,  in  conditions  which  are  either  exceptionally  good  or  axception- 

(4) 

ally  bad,  all  thrae  methods  fail  in  reliable  range  prediction.  The 
implicit-explicit  method  has  a  potential  which  is  not  available  with 
current  methods  for  range  prediction  under  the  conditions  of  a  variable 


local  environment. 
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Appendix 

COMPUTER  PROGRAM 

S0UNDW  is  a  main  routine  written  in  FORTRAN  IV  for  use  on  the  IBM- 
7044,  but  it  also  may  be  used  on  the  various  machines  capable  of  using 
the  FORTRAN  language . 

Input,  read  in  by  FORMAT  (15),  consists  f: 

1.  Length  of  grid  -  MDIM 

2.  Width  of  grid  -  NMAX 

3 .  Expanding  boundary  for  M  -  MMAX 

4.  Maximum  number  of  time  steps  -  MAXST 

5.  Initial  time  step  number  -  NST 

6.  Frequency  -  F 

7.  Length  of  grid  interval  -  AL 

8.  Velocity  -  VEL 

9.  Time  step  size  -  AT 

10.  Rho  -  RH0 

11.  Gamma  -  GAM 'A 

12.  Amplitude  -  AMPU 

Output  is  stored  on  F0RTPAN  1/0  unit  9 

1.  SEP,  SEPP  -  density 

2.  UP,  UPP  -  X-component  of  velocity 

3.  WP,  WPP  -  Z-component  of  velocity 

where  P  refers  to  the  first  half-time  step,  and  PP  refers  to  the  second 
half-time  step. 


on  on  ooonooononoonnw 
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8FTC  SOUNOW  REF 

SIMULATES  SOUND  PROPAGATION  UNDERWATER 

********** 

♦  SOUNCW  * 

********** 

M.C  •  SMIiH 
THE  RANO  CORP. 

SANTA  MON ICA, CAL  IF • 

WRITTEN  BY  R.H.  MAYALl,  JULY, 1967 
REVISEO  BY  S.V.  HUBER.  TECEMBER, 1967 

DIMENSION  PO IT ( 20,  124 1 

D  IMF  *1S  ION  U(20,2?2),SE(20,222),W(20*  222 ), WPI 20,22? ) 
DIMENSION  SEP! 20, 222), SEPPI 20,222! , UPI 20,222) ,UPPI 20,222) , 
1  WPP ( 20, 222 ) 

EQUIVALENCE  ( SE ( 1,  1 ) , SEP( l, 1 ) , SEPP 1 1 , 1 ) ) 

EQUIVALENCE  (U( 1, 1 ) ,UP( 1, 1 ) ) 

EQUIVALENCE  ( W (  1, 1 ) , UPP ( 1  *  1 ) ) 

EQUIVALENCE  (WPI  1, 1) ,WPP( 1,1)) 

DIMENSION  A(222),B(222),P(222),Q(222),R(222),S(222) 
DIMENSION  aOAV(222), X1DAV(  222) 

DIMENSION  CSQ121I 
DIMENSION  DES(6)  .DESK  2) 

DIMENSION  N0( 222 ) 

OATA  OES  / 3H  U, 4H  UP, AH  SE,AH  SEP,3H  W,AH  WP/ 

DATA  0ES1  /5HF IRST, 6HSEC0ND/ 

READ ( 5, AOA )  MOIM 
READ ( 5, AOA  )  NM AX 
REAO (5, AOA )  MMAX 
READ ( 5, AOA )  MAXST 
READI5, AOA)  NST 
READ ( 5, A06 )  F 
REAO (5, 406)  AL 
REAO (5, A06  )  VCL 
READ ( 5, A06 )  AT 
READ! 5, A06 )  RHO 
READ ( 5, A06  )  GAMMA 
READ (5,406)  AMPU 
AOA  FORMAT ( 15  ) 

A06  FORMAT (F10. A) 

DO  10  I=1,MDIM 
10  NOCI)  *  I 

*♦  MMAX  SHOULD  BE  EVEN  INITIALLY. 

NMAX1  *  NMAX-1 
MMAX 1  *  MMAX- 1 
AL  *  AL*30.5 
CF  *  .975914506 
VEL  *  V EL* 30.5 
Cl  «  AT/AL 
C2  *  Cl/AL 

INITIAL  CONDITIONS 
DO  20  M» 1,M0IM 
AIM)  *  0. 

B(M)  *  0. 


o  o  no  no  no  no  no 
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P(M)  *  0. 

0(M)  =  0. 

RIM  )  x  0. 

S(  M )  *  0. 

00  20  N*  1 1  NM'.X 
U(NtM)  *  0. 

SEtNtM)  *  0.0 
W  (N»  M  )  -  G. 

20  WP(N,M.  *  0. 

PI  *  3. 14 15926 
FPI  «  4.0*PI*F*AT 
TR0L  *  CI*RH0 
PL2  *  CI/AL/RH0 
PLGAM  *  PL2*GAMMA 
00  30  N*2t  NMAX 
CSQ(N)  *  Cl*V£L**2/RH0 
30  CONTINUE 

XDAVf XIDAV  USED  ONLY  FOR  A  CYLINDRICAL  COORDINATE  SYSTEM 
DO  40  I*2» MDIM 
XM  *  AL*FL0AT 1  I- I ) 

XMM*  AL*FL0AT< 1-2) 

XDAV(I)  *  2.*XM/( XM*XMM )*TRCL 
XIDAV(I)  *  2 .*XMM/ ( XM+XMM ) ♦TROL 
40  CONTINUE 

FIRST  HALF-TIME  STEP  CALCULATIONS  BEGIN 
50  CONTINUE 

EXPLICIT  CALCULATION  OF  WP 

XNST  *  NST 

DO  70  M*2»MMAX 

MM  *  M-l 

MMM  «  M*1 

DO  fcO  N=2» NMAXl 

NN  *  N-l 

NNN  =  N*l 

HOLDI  *  W ( N»  M )  -  C  SQ ( NNN )  *  aE  ( NNN  *  M )  ♦  C  SQ( N) *SE ( N«  M) 
WP!N*M)  *  HOLDI  ♦  PLGAM*! W(NNNyM)  -  2.0*W(N,M)  ♦  W(NN,M)) 
60  CONTINUE 

WP(NMAXyM)  *  0.0 
70  CONTINUE 

RIN*  L ) « S( N» 1 )  DEFINE  BURY.  CONDS.  FOR  UP 
RID  =  0.0 
L  *  MMAX 

IMPLICIT  CALCULATION  OF  P,Q,RfS 
DO  100  N*  2»NMAX 
S(l)  «  0.0 

POINT  SOURCE 

IFIN.EQ.14)  Sill  *  AMPU*SIN! FPI* XNST) 

NN  *  N-l 
00  80  M*2,L 
MM  =  M-l 
MMM  *  M*l 

A(M )  «  SE(NfM)  ♦  TROL* ( W( NN» M )  -  W(N,M) ) 

POEM  «  1.0  ♦  TROL*R(MM) 

P(M)  *  TROL/PDEM 


o  o  o  o  no  on  on  no 
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0(H)  *  I  A(  M )  ♦  TR0L*S(MM) )/PDEM 

IFIM. EQ.MMAX)  GO  TO  80 

BIM  )  *  U( N,M  )  ♦  PLGAM*(  U(  N,  ’IMM)  -2.0*UlN,M)  ♦  U( N, MM)  ) 

RNUH  >  CSQIN  ) 

RDEM  «  1.0  ♦  RNUM*P( M ) 

RIM)  >  RNUM/RDEM 
SIM)  *  ( B ( M  )  ♦  RNUM*Q(M) ) /RDEM 
80  CONTINUE 

UP(N,L )*  0.0 

M*L  - 

IMPLICIT  CALCULATION  OF  SEP, UP 
DO  90  J* 2»MMAX 
MM  «  M-l 

SEP (N»M )  »-P(M )*UP(N,M)+Q(M) 

UP( N, MM )  «-K(MM)*SEP(N,M)+S(MM) 

90  M  «  M-l 
100  CONTINUE 
MX  *  MMAX 
ASSIGN  110  TO  KK 
GO  TO  190 

SECONO  HALF-TIMF  STEP  CALCULATIONS  BEGIN 
110  CONTINUE 

EXPLICIT  CALCULATION  OF  UPP 

DO  130  N*2,NMAX 

NN  «  N-l 

NNN  *  N*l 

DO  120  M«2,MMAXl 

MM  -  M-l 

MMM  *  M*1 

UPP(N,M)*UP(N,M)  -  CSQ(N)*(SEP(N,MMM)  -  SEP(N,M)|  ♦  PLGAM* 

1  (UP  ( N,MMM  )— 2.0*UP(N,  M I  ♦  UPIN,MN)) 

120  CONTINUE 
130  CONTINUE 

BOUNOARY  CONDITIONS  ON  UPP 
DO  l AO  N*-*  1,NMAX 
UPP IN* 1 )  *  0.0 
140  UPP I N, MM AX )  x  0.0 

UPP I  14,1)  ■  AMPU*SIN(FPI4(XNST4.5) ) 

IMPLICIT  CALCULATION  OF  P,Q,R,S 
DO  170  M* 2, MMAX 

R(1,M),S( 1 , M )  DEFINE  BDRY.  CONDS.  FOR  WPP 

R(  1  )  *  0.0 

SI  l  )  *  o.O 

MM  «  M-l 

MMM  «  M»1 

DO  150  N*2,NMAX 

NN  «  N-l 

NNN  *  N+l 

AIN)  »  SEP(N,M)  -  TROL* I  UP  I N ,M )-UP( N ,MM ) ) 

POE  «  1.0  ♦TRCL*R I NN  ) 

PIN)  xTROL/PDE 
0  IN )  « I A ( N )  ♦TROL*S(NN) )/PDE 

IFIN. EQ.NMAX)  GO  TO  150 

BIN )*WP (N, M ) ♦  PLGAM*(WP(NNN,M)-2.*WP(N,M)+WP(NN,M) ) 


o  o  o  o 
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RDE  *  1.0  ♦  CSC,!N)*P!N) 

R IN )  *  CSQ!NNN)/RDE 
S(N)  *  I  B!  N )  ♦  Q!N)*CSQIN) )/RDE 
150  CONTINUE 

CONI  *  CSQiNMAX) 

N*  NMAX 
C 

C  BOUNDARY  CONDITIONS  ON  WPP 
WPP ( NMAX, M )  *  0.0 

IMPLICIT  CALCULATION  OF  SEPP.WPP 
DO  160  J»2,NMAX 
NN  *  N-l 

SEPP(NfM)  *  -P(N )*WPP! N»M)  ♦  Q(N) 

WPP( NN» M I  *  -R!NN)*SEPPIN,M)  4  SINN) 

160  N  *  N-l 
170  CONTINUE 

DO  180  J* 1*  MMAX 
DO  180  I* I* NMAX 
U(ItJ)  *  UPPII.J) 

180  W!I,J)  *  WPPII.J) 

ASSIGN  200  TO  KK 
190  CONTINUE 

DO  2  N*  It  20 
DO  3  M*  1,124 

IF(SE(N,M).GT.POITIN,M))  P0IT1N,M)  *  SEIN.M) 
3  CONTINUE 
2  CONTINUE 

GO  TO  KK, ( 110*200) 

200  CONTINUE 

NST  *  NST ♦! 


EXPANDING  BOUNDARY 
IF(MMAX.LT.MCIM)  MMAX  *  MMAX+2 
MMAX 1  *  MMAX-l 
IFINST.GT.MAXST)  GO  TO  210 
GO  TO  50 
210  CONTINUE 

WRITE19)  POIT 
DO  13  M*  1,126 

PRINT  16,  ( PO I T ( N, M ) ,  N«l,20) 

13  CONTINUE 
STOP 

16  FORMAT  (1H  , 10E13.5) 

300  FORMAT  1 1H  , 13, ?X , 1 IE  1 1.6 ) 

301  FORMAT! 1H1,6X,6HTHE  ,A6,9HVALUES  AT,1X,A6,17H  HALF  OF  TIMESTFP, 15  ) 
600  FORMAT! IH  ,3X, 13) 

602  FORMAT! 1H0,6HNST  =  , I6.23HALL  VALUES  ARE  AVERAGED) 

END 


♦ENTRY 

222 

20 

60 

350 

1 

20.0 

12.5 

5000.0 

.00125 

1.025 


SOUNDW 

LENGTH  OF  GRID 
WIDTH  OF  GRID 
EXPANDING  BOUNDRY  FOR  M 
MAX.  NO.  OF  TIMESTEPS 
INITIAL  TIMESTEP  NO. 
FREQUENCY 

LENGTH  OF  GRID  INTERVAL 

VELOCITY 

TIMESTEP  SIZE 

RHO 


0114 

6.0 


GAMMA 

AMPLITUDE 
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SIBSYS 

SCLOSE  S.SU06, REMOVE 

$ IBFTC  PLOTP 

DIMENSION  AMQDESI 200) 

DIMENSION  P0ITI20, 124) 

DIMENSION  SIZE(4)»X(2500)»Y(2500) 

DATA  SIZE/. 75*1. 0* 1.25(1.5/ 

READ  (9)  POIT 

CALL  MODE SGI AMODESf 0) 

CALL  SETSMG  ( AMODESt 84,6HXXXXXX) 

CALL  SETSMG  ( AMODESt 94( 1. ) 

CALL  OBJCTGI AMODESf 0.t0.f 1.333333* .3) 

CALL  SU8JEG  (AMODESt 0.0(0.0*124. ,20. ) 

CALL  GRIDG  I  ANODES, 1.0, 1.0,0, 0) 

DO  20  1*1,4 
IX  *  0 

DO  20  K*l, 20 

DO  40  L*l, 124 

GO  TC  (60,70,80,90),  I 

60  IE(POIT(K,L).LT.  . 50E-5.AND. POI T( K,L ) .GE.  .40E-5)  GO  TO  50 
GO  TO  40 

70  IF(POIT(K,L).LT.  .40E-5.AND.P0I T( K,L) .GE.  .30E-5)  GO  TO  50 
GO  TO  40 

80  IF(POIT(K,L).LT.  . 30E-5. AND. POI T(K,L).Gt.  .20E-5)  GO  TO  50 
GO  TO  40 

90  IF(POIT(K,L).LT.  . 20E-5. AND.POI T(K,L).GE.  .OOE-1)  GO  TO  50 
GO  TO  40 
50  IX  *  IX*1 
X(rX)  *  L 
Y(IX)  »  K 
40  CONTINUE 
30  CONTINUE 

CALL  SETSMG( ANODES, 53, SIZE ( I ) ) 

CALL  POINTG  (ANODES, IX, X, Y) 

20  CONTINUE 

CALL  PICTRGI ANODES, 1,0,1) 

CALL  EXITG  (ANODES) 

CALL  EXIT 
END 
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